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ABSTRACT 

Pre-main-sequence stars are observed to be surrounded by both accretion flows and some kind of wind or 
jet-like outflow. Recent work by Matt and Pudritz has suggested that if classical T Tauri stars exhibit stellar 
winds with mass loss rates about 0.1 times their accretion rates, the wind can carry away enough angular mo- 
mentum to keep the stars from being spun up unrealistically by accretion. This paper presents a preliminary 
set of theoretical models of accretion-driven winds from the polar regions of T Tauri stars. These models are 
based on recently published self-consistent simulations of the Sun's coronal heating and wind acceleration. In 
addition to the convection-driven MHD turbulence (which dominates in the solar case), we add another source 
of wave energy at the photosphere that is driven by the impact of plasma in neighboring flux tubes undergo- 
ing magnetospheric accretion. This added energy, determined quantitatively from the far-field theory of MHD 
wave generation, is sufficient to produce T Tauri-like mass loss rates of at least 0.01 times the accretion rate. 
While still about an order of magnitude below the level required for efficient angular momentum removal, 
these are the first self-consistent models of T Tauri winds that agree reasonably well with a range of obser- 
vational mass loss constraints. The youngest modeled stellar winds are supported by Alfven wave pressure, 
they have low temperatures ("extended chromospheres"), and they are likely to be unstable to the formation of 
counterpropagating shocks and clumps far from the star. 

Subject headings: accretion, accretion disks — stars: coronae — stars: mass loss — stars: pre-main sequence 
— stars: winds, outflows — turbulence 



1. INTRODUCTION 

Our current state of knowledge about how stars and plan- 
ets are formed comes from an intertwined web of observa- 
tions (spanning the electromagnetic spectrum) and theoretical 
work. The early stages of low-mass star formation comprise 
a wide array of inferred physical processes, including disk ac- 
cretion, various kinds of outflow, and magnetohydrodynamic 
(MHD) activity on time scales ranging from hours to millen- 
nia (see, e.g., Lada 1985; Bertout 1989; Appenzeller & Mundt 
1989; Hartmann 2000; Konigl & Pudritz 2000; McKee & Os- 
triker 2007; Shu et al. 2007). A key recurring issue is that 
there is a great deal of mutual interaction and feedback be- 
tween the star and its circumstellar environment. This inter- 
action can determine how rapidly the star rotates, how active 
the star appears from radio to X-ray wavelengths, and how 
much mass and energy the star releases into its interplanetary 
medium. 

A key example of circumstellar feedback is the mag- 
netospheric accretion paradigm for classical T Tauri stars 
(Lynden-Bell & Pringle 1974; Uchida & Shibata 1984; Ca- 
menzind 1990; Konigl 1991). Because of strong (~1 kG) 
stellar magnetic fields, the evolving equatorial accretion disk 
does not penetrate to the stellar surface, but instead is stopped 
by the stellar magnetosphere. Accretion is thought to proceed 
via ballistic infall along magnetic flux tubes threading the in- 
ner disk, leading to shocks and hot spots on the surface. The 
primordial accretion disk is dissipated gradually as the star en- 
ters the weak-lined T Tauri star phase, with a likely transition 
to a protoplanetary dust/debris disk. 

Throughout these stages, solar-mass stars are inferred to ex- 
hibit some kind of wind or jet-like outflow. There are several 
possible explanations of how and where the outflows arise, in- 
cluding extended disk winds, X-winds, impulsive (plasmoid- 
like) ejections, and "true" stellar winds (e.g., Paatz & Camen- 



zind 1996; Calvet 1997; Konigl & Pudritz 2000; Dupree et 
al. 2005; Edwards et al. 2006; Ferreira et al. 2006; Gomez de 
Castro & Verdugo 2007; Cai et al. 2008). Whatever their ori- 
gin, the outflows produce observational diagnostics that indi- 
cate mass loss rates exceeding the Sun's present mass loss rate 
by factors of 10 3 to 10 6 . It is of some interest to evaluate how 
much of the observed outflow can be explained solely with 
stellar winds, since these flows are locked to the star and thus 
are capable of removing angular momentum from the system. 
Recent work by Matt & Pudritz (2005, 2007, 2008) has sug- 
gested that if there is a stellar wind with a sustained mass loss 
rate about 10% of the accretion rate, the wind can carry away 
enough angular momentum to keep T Tauri stars from being 
spun up unrealistically by the accretion. 

Despite many years of study, the dominant physical pro- 
cesses that accelerate winds from cool stars have not yet been 
identified conclusively. For many stars, the large-scale en- 
ergetics of the system — i.e., the luminosity and the gravita- 
tional potential — seem to determine the overall magnitude of 
the mass loss (Reimers 1975; 1977; Schroder & Cuntz 2005). 
Indeed, for the most luminous cool stars, radiation pressure 
seems to provide a direct causal link between the luminos- 
ity L* and the mass loss rate M w j n d (e.g., Gail & Sedlmayr 
1987; Hofner 2005). However, for young solar-mass stars, 
other mediating processes (such as coronal heating, waves, or 
time-variable magnetic ejections) are more likely to connect 
the properties of the stellar interior to the surrounding out- 
flowing wind. For example, magnetohydrodynamic (MHD) 
waves have been studied for several decades as a likely way 
for energy to be transferred from late-type stars to their winds 
(Hartmann & MacGregor 1980; DeCampli 1981; Airapetian 
et al. 2000; Falceta-Goncalves et al. 2006; Suzuki 2007). 

In parallel with the above work in improving our under- 
standing of stellar winds, there has been a great deal of 
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FIG. 1. — Summary sketch of the accretion and wind geometry discussed 
in this paper. (1) Magnetospheric accretion streams are assumed to impact 
the star at mid-latitudes. (2) Dense clumps in the accretion streams generate 
MHD waves that propagate horizontally over the stellar surface. (3) Polar 
magnetic field lines exhibit enhanced photospheric wave activity and thus 
experience stronger coronal heating and stellar wind acceleration. 

progress toward identifying and characterizing the processes 
that produce the Sun 's corona and wind. It seems increasingly 
clear that closed magnetic loops in the low solar corona are 
heated by small-scale, intermittent magnetic reconnection that 
is driven by the continual stressing of their footpoints by con- 
vective motions (e.g., Aschwanden 2006; Klimchuk 2006). 
The open field lines that connect the Sun to interplanetary 
space, though, appear to be energized by the dissipation of 
waves and turbulent motions that originate at the stellar sur- 
face (Tu & Marsch 1995; Cranmer 2002; Suzuki 2006; Kohl 
et al. 2006). Parker's (1958) classic paradigm of solar wind 
acceleration via gas pressure in a hot (T ~ 10 6 K) corona 
still seems to be the dominant mechanism, though waves and 
turbulence have an impact as well. A recent self-consistent 
model of turbulence-driven coronal heating and solar wind 
acceleration has succeeded in reproducing a wide range of 
observations with no ad-hoc free parameters (Cranmer et al. 
2007). This progress on the solar front is a fruitful jumping- 
off point for a better understanding of the basic physics of 
winds and accretion in young stars. 

The remainder of this paper is organized as follows. § 2 
presents an overview of the general scenario of accretion- 
driven MHD waves that is proposed here to be important for 
driving T Tauri mass loss. In § 3 the detailed properties of an 
evolving solar-mass star are presented, including the funda- 
mental stellar parameters, the accretion rate and disk geom- 
etry, and the properties of the clumped gas in the magneto- 
spheric accretion streams. These clumped streams impact the 
stellar surface and create MHD waves that propagate horizon- 
tally to the launching points of stellar wind streams. § 4 de- 
scribes how self-consistent models of these wind regions are 
implemented, and § 5 gives the results. Finally, § 6 contains 
a summary of the major results of this paper and a discus- 
sion of the implications these results may have on our wider 
understanding of low-mass star formation. 

2. BASAL WAVE DRIVING FROM INHOMOGENEOUS ACCRETION 

Figure 1 illustrates the proposed connection between ac- 
cretion and stellar wind driving that is explored below. The 
models make use of the idea that magnetospheric accretion 
streams are likely to be highly unstable and time-variable, 
and thus much of the material deposited onto the star is ex- 
pected to be in the form of intermittent dense clumps. The 
impact of each clump is likely to generate MHD waves on 
the stellar surface (Scheurwater & Kuijpers 1988) and these 
waves can propagate across the surface, spread out the ki- 
netic energy of accretion, and inject some of it into the global 



magnetic field. The models simulate the time-steady prop- 
erties of the polar magnetic flux tubes, wherein the source 
of wave/turbulence energy at the photosphere comes from 
two sources: (1) accretion-driven waves that have propagated 
over the surface from the "hot spots," and (2) the ever-present 
sub-photospheric convection, which shakes the magnetic field 
lines even without accretion, and is the dominant source of 
waves for the present-day Sun. 

There is substantial observational evidence that the accre- 
tion streams of classical T Tauri stars are clumped and in- 
homogeneous (see, e.g., Gullbring et al. 1996; Safier 1998; 
Bouvier et al. 2003, 2004, 2007; Stempels 2003). Strong vari- 
ations in the accretion signatures take place over time scales 
ranging between a few hours (i.e., the free-fall time from the 
inner edge of the disk to the star) and a few days (i.e., the 
stellar rotation period). Some of this variability may arise 
from instabilities in the accretion shock (Chevalier & Ima- 
mura 1982; Mignone 2005), with the possibility of transitions 
between stable and unstable periods of accretion (e.g., Ro- 
manova et al. 2008). Observations of the bases of accretion 
streams indicate the presence of turbulence (Johns & Basri 
1995) as well as discrete events that could signal the presence 
of magnetic reconnection (Giardino et al. 2007). Alternately, 
some of the observed variability may come from larger-scale 
instabilities in the torqued magnetic field that is threaded by 
the disk (e.g., Long et al. 2007; Kulkarni & Romanova 2008) 
and which could also drive periodic reconnection events (van 
Ballegooijen 1994). The latter variations would not necessar- 
ily be limited to the rotation time scale, since "instantaneous" 
topology changes can happen much more rapidly. 

The mechanisms by which impulsive accretion events (i.e., 
impacts of dense clumps) can create MHD waves on the stel- 
lar surface are described in more detail in § 3.4. This paper 
uses the analytic results of Scheurwater & Kuijpers (1988) to 
set the energy flux in these waves, based on the available ki- 
netic energy in the impacts. It is useful to point out, however, 
that there are relatively well-observed examples of impulse- 
generated waves on the Sun that help to justify the overall 
plausibility of this scenario. These are the so-called "More- 
ton waves" and "EIT waves" that are driven by solar flares 
and coronal mass ejections (CMEs). Moreton waves are lo- 
calized ripples seen in the Ha chromosphere expanding out 
from strong flares and erupting filaments at speeds of 500 to 
2000 km s -1 (Moreton & Ramsey 1960). These seem to be 
consistent with fast-mode MHD shocks associated with the 
flaring regions (e.g., Uchida 1968; Balasubramaniam et al. 
2007). A separate phenomenon, discovered more recently in 
the extreme ultraviolet and called EIT waves (Thompson et 
al. 1998), is still not yet well understood. These waves prop- 
agate more slowly than Moreton waves (typically 100 to 400 
km s" 1 ), but they are often seen to traverse the entire diam- 
eter of the Sun and remain coherent. Explanations include 
fast-mode MHD waves (Wang 2000), solitons (Wills-Davey 
et al. 2007), and sheared current layers containing enhanced 
Joule heating (Delannee et al. 2008). These serve as ample 
evidence that impulsive phenomena can generate MHD fluc- 
tuations that travel across a stellar surface. 

Finally, it is important to clarify some of the limitations of 
the modeling approach used in this paper. The models include 

' A somewhat similar wave-generation scenario was suggested by Vascon- 
celos et al. (2002) and Elfimov et al. (2004), but their model was applied only 
to the energization of the accretion streams themselves, not to the regions 
exterior to the streams. 
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only a description of the magnetospheric accretion streams 
(with an assumed dipole geometry) and the open flux tubes 
at the north and south poles of the star. Thus, there is no at- 
tempt to model either disk winds or the closed-field parts of 
the stellar corona, despite the fact that these are likely to con- 
tribute to many key observational diagnostics of T Tauri cir- 
cumstellar gas. In addition, the accretion streams themselves 
are included only for their net dynamical impact on the polar 
(non-accreting) regions, and thus there is no need to describe, 
e.g., the temperature or ionization state inside the streams. 
The wind acceleration in the polar flux tubes is treated with a 
one-fluid, time-steady approximation similar to that described 
by Cranmer et al. (2007) for the solar wind. To summarize, 
this paper is a "proof of concept" study to evaluate how much 
about T Tauri outflows can be explained with only polar stellar 
winds that are energized by accretion-driven waves. 

3. BASELINE SOLAR-MASS EVOLUTION MODEL 
3.1. Stellar Parameters and Accretion Rate 

To begin exploring how accretion-driven waves affect the 
winds of young stars, a representative evolutionary sequence 
of fundamental parameters was constructed for a solar-type 
star. The adopted parameters are not intended to reproduce 
any specific observation in detail and should be considered il- 
lustrative. The STARS evolution code 2 (Eggleton 1971, 1972, 
1973; Eggleton et al. 1973; Han et al. 1994; Pols et al. 1995) 
was used to evolve a star having a constant mass M* = 1 M Q 
from the Hayashi track to well past the current age of the Sun. 
Neither mass accretion nor mass loss were included in the 
evolutionary calculation, and a standard solar abundance mix- 
ture (X = 0.70, Z = 0.02) was assumed. All other adjustable 
parameters of the code were held fixed at their default values. 

Figure 2 presents a summary of the modeled stellar param- 
eters. A Hertzsprung-Russell (H-R) diagram is shown in Fig- 
ure 2a, with the bolometric luminosity L* plotted as a func- 
tion of the effective temperature T e ff for both the evolutionary 
model (with a subset of stellar ages t indicated by symbols) 
and an approximate main sequence band for luminosity class 
V stars (de Jager & Nieuwenhuijzen 1987). The current age 
of the Sun is denoted by t = 4.6 x 10 9 yr, or log I0 f = 9.66. 
Figure 2b shows the age dependence of a selection of relevant 
length scales, including the stellar radius and the photo- 
spheric pressure scale height //» . The latter quantity is defined 
assuming a pure hydrogen gas, as 



mu GM* 



(1) 



where is Boltzmann's constant, »?h is the mass of a hydro- 
gen atom, and G is the Newtonian gravitation constant. Note 
that for most of the applications below, only the relative vari- 
ation of with age is needed and not its absolute value. The 
scale height is assumed to be proportional to the horizontal 
granulation scale length (e.g., Robinson et al. 2004), which in 
turn governs the perpendicular correlation length of Alfvenic 
turbulence in the photosphere (see § 4). 

The other major parameter to be specified as a function of 
age is the mass accretion rate M acc . For classical T Tauri stars, 
it is generally accepted that the accretion rate decreases with 
increasing stellar age t , but there is a large spread in measured 
values that may be affected by observational uncertainties in 
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FIG. 2. — Basic properties of 1 Mq star, (a) H-R diagram showing the main 
sequence (gray band) and the modeled evolutionary track, with several rep- 
resentative ages (open circles) including the present-day Sun (filled circle), 
(b) Length scales, plotted in units of solar radii, as a function of age: stellar 
radius (solid line), inner edge of accretion disk (dashed line), accretion clump 
radius (dot-dashed line), and photospheric pressure scale height (dotted line). 

both M acc and t. In order to determine a representative age 
dependence for M acc (f), we utilized tabulated accretion rates 
from Hartigan et al. (1995) and Hartmann et al. (1998), who 
interpreted optical/UV continuum excesses as an "accretion 
luminosity" that comes from the kinetic energy of accreted 
gas impacting the star. There is some overlap in the star lists 
from these two sources, and differences in the diagnostic tech- 
niques resulted in different values of both M acc and t for some 
stars common to both lists. Both values have been retained 
here, and thus Figure 3a shows a total of 88 data points from 
both lists. 3 The observational uncertainties are not shown; 
they may be as large as an order of magnitude in M acc and a 
factor of ^3 in the age (see, e.g., Figure 3 of Muzerolle et al. 
2000). 

Highlighted in Figure 3a, as filled symbols, are a subset of 
stars that appear to be on the 1 M Q Hayashi track (i.e., they 
have T e ff between about 4000 and 4500 K). An initial attempt 
to fit these 35 stars to a power-law age dependence yielded 



M acc w 2.8 x lO^Moyr -1 



10 6 yr 



-1.22 



(2) 



2 The December 2003 
http://www.ast.cam.ac.uk/~stars/ 



version was obtained from: 



3 This figure shows relatively nearby Galactic stars only. In the LMC, the 
mass accretion rates seem to be much higher at ages around 10 Myr (e.g., 
Romaniello et al. 2004), which could indicate a substantial metallicity de- 
pendence for many properties of the T Tauri phase. 
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FIG . 3 . — Mass accretion rates and mass loss rates for classical T Tauri stars. 
(a) Measured accretion rates from Hartigan et al. (1995) and Hartmann et al. 
(1998) (open circles) with the subset of roughly solar-mass stars highlighted 
(filled circles). The best-fit curve from eq. {5J is also shown (solid line) as 
well as two measured rates for the stars TW Hya and MP Mus (error bars); 
see text for details, (b) Ratio of mass loss rate to mass accretion rate for 
measured T Tauri stars, where the symbols have the same meanings as in 
panel (a). 

However, the theoretical accretion disk models of Hartmann 
et al. (1998) found that the exponent rj in M acc oc r n is most 
likely to range between about 1.5 and 2.8, with the lower end 
of that range being most likely. Thus, the fitted value of 77 = 
1.22 above was judged to be too low, and a more realistic fit 
was performed by fixing 77 = 1 .5 and finding 



M acc w 2.2 x 10" s M o yr" 



10 6 yr 



-1.5 



(3) 



which is used in the wind models below and is plotted in Fig- 
ure 3a. The scatter around the mean fit curve has a roughly 
normal distribution (in the logarithm of M acc ) with a standard 
deviation ranging between a factor of 6 (for the 35 solar-mass 
stars) to 8 (for all 88 data points) on either side of the repre- 
sentative trend given by equation ||3}. The sensitivity of the 
resulting accretion-driven wind to this scatter is explored fur- 
ther below in § 5.4. 

Classical T Tauri magnetospheric accretion is believed to 
end at an age around 10 Myr (e.g., Strom et al. 1989; Bou- 
vier et al. 1997). This is probably coincident with the time 
when the inner edge of the accretion disk grows to a larger 
extent than the Keplerian corotation radius, with "propeller- 
like" outflow replacing the accretion (Illarionov & Sunyaev 



1975; Ustyugova et al. 2006). Rotation is not included in the 
models presented here, so this criterion is not applied directly. 
However, for advanced ages (e.g., t > 100 Myr) equation (01 
gives increasingly weak accretion rates that end up not having 
any significant effect on the stellar wind. Thus, even though 
the Figures below plot the various accretion-driven quantities 
up to an age of 1 Gyr (log 10 f = 9), the models below do not 
apply any abrupt cutoff to the accretion rate. 

Figure 3b plots mass loss rates M w i n d for a subset of classi- 
cal T Tauri stars that have measurements of blueshifted emis- 
sion in forbidden lines such as [O I] A6300 (Hartigan et al. 
1995). The mass loss rates are shown as an efficiency ratio 
M wl!U i/M acc which depends on the combined uncertainties of 
both the outflow and accretion properties. The largest ratio 
shown, at log 10 f « 6.9, is for the jet of HN Tau, and the orig- 
inal value from Hartigan et al. (1995) has been replaced with 
the slightly lower revised value from Hartigan et al. (2004). It 
is uncertain whether the measured outflows originate on the 
stellar surface or in the accretion disk (see, e.g., Paatz & Ca- 
menzind 1996; Calvet 1997; Ferreira et al. 2006), but these 
rates can be used as upper limits for any stellar wind compo- 
nent. 

Figure 3 also gives additional information about two of the 
oldest reported classical T Tauri stars: TW Hya and MP Mus. 
Note, however, that their extremely uncertain accretion rates 
were not included in the fitting for M acc (f). The age of TW 
Hya is quite uncertain; it is often given as ^10 Myr (e.g., 
Muzerolle et al. 2000), but values as large as 30 Myr have 
been computed (Batalha et al. 2002) and the plotted error 
bar attempts to span this range. Similarly, the age of MP 
Mus seems to be between 7 and 23 Myr (e.g., Argiroffi et 
al. 2007). The plotted upper and lower limits on M acc for TW 
Hya come from from Batalha et al. (2002) and Muzerolle et 
al. (2000), respectively, and the range of values for M w i„d in 
Figure 3b were taken from Dupree et al. (2005). The plotted 
ratio, however, divides the observational limits on M w j nl j by 
a mean value of M acc = 1 .4 x 10~ 9 M Q yr" 1 , in order to avoid 
creating an unrealistically huge range. For MP Mus, the mean 
accretion rate of 5 x 10 -11 M Q yr" 1 was derived from X-ray 
measurements by Argiroffi et al. (2007). The plotted error 
bars forM acc were estimated by using the la uncertainties on 
the X-ray emission measure and column density. 



3.2. Magnetospheric Accretion Streams 

The dynamical properties of accretion streams are modeled 
here using an axisymmetric dipole magnetic field, coupled 
with the assumption of ballistic infall (e.g., Calvet & Hart- 
mann 1992; Muzerolle et al. 2001). Although it is almost 
certain that the actual magnetic fields of T Tauri stars are not 
dipolar (Donati et al. 2007; Jardine et al. 2008), this assump- 
tion allows the properties of the accretion streams to be calcu- 
lated simply and straightforwardly as a function of evolution- 
ary age. 

Stellar field lines that thread the accretion disk are bounded 
between inner and outer radii r m and r out as measured in the 
equatorial plane. We assume the inner radius — also called the 
"truncation radius" — is described by Konigl's (1991) applica- 
tion of neutron-star accretion theory (see also Davidson & Os- 
triker 1973; Eisner & Lamb 1977; Ghosh & Lamb 1979a,b). 
This expression for r m is given by determining where the mag- 
netic pressure of the inner dipole region balances the gas pres- 
sure in the outer accretion region. Assuming free-fall in the 
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accretion stream, 



r m = Pcl 



( BtRl 



1/7 



(4) 



where the scaling constant /?gl describes the departure from 
ideal magnetostatic balance; i.e., /3gl = 1 gives the standard 
"Alfven radius" at which the pressures balance. Following 
Konigl (1991), the value /3gl = 0.5 is used here. The outer 
disk radius r out may be as large as the Keplerian corotation 
radius, but it may not extend that far in reality (e.g., Hartmann 
et al. 1994; Bessolaz et al. 2008). For simplicity, the outer 
disk radius is assumed to scale with the inner disk radius as 



r in (l+e) 



(5) 



where a constant value of e = 0. 1 was adopted. This value falls 
within the range of empirically constrained outer/inner disk 
ratios used by Muzerolle et al. (2001) to model Ha profiles; 
they used effective values of e between 0.034 and 0.36. The 
intermediate value of 0. 1 produces reasonable magnitudes for 
the filling factors of accretion stream footpoints on the stellar 
surface (see below). 

Note that equation (0]i above requires the specification of 
the surface magnetic field strength B*. In the models of the 
open flux tubes (in the polar regions) used below, a solar- 
type magnetic field is adopted. The photospheric value of 
fi» pa 1500 G is held fixed for the footpoints of the stellar wind 
flux tubes (see Cranmer & van Ballegooijen 2005). For the 
mid-latitude field at the footpoints of the accretion streams, 
though, a slightly weaker field of 1000 G is used. 

Figure 2b shows the resulting age dependence for r; n . For 
the youngest modeled "protostars" (t < 10 4 yr), the accretion 
rate is so large that r m would be smaller than the stellar radius 
itself. In that case, the accretion disk would penetrate down to 
the star and there would be no magnetospheric infall. Thus, 
the youngest age considered for the remainder of this paper is 
t = 13.5 kyr (i.e., log 10 f =4.13), for which r; n = R*. Over most 
of the classical T Tauri age range (0.1-10 Myr), r m decreases 
slightly in absolute extent, but the ratio r m /R* increases and 
is approximately proportional to ?° 2 . 

For any radius r > R* in the equatorial disk, an aligned 
dipole field line impacts the stellar surface at colatitude 9, 
where R* = rsin 2 #. This allows the colatitudes 6 ln and 9 oM 
to be computed from r m and r out , and these can be used to 
compute the fraction of stellar surface area S that is filled by 
accretion streams, with 



S = COS# nu t-COS#i, 



(6) 



(see also Lamzin 1999). Note that both the northern and 
southern polar "rings" are taken into account in the above ex- 
pression. The filling factor 6 depends sensitively on the geom- 
etry of the accretion volume; i.e., for a dipole field, it depends 
only on the outer/inner radius parameter e and the ratio r m /R*. 
Similarly, the fractional area subtended by both polar caps, to 
the north and south of the accretion rings, is given by 

l-cos0 out . (7) 

and estimated it to 
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Dupree et al. (2005) called this quantity 
be -0.3 for TW Hya. 

Accreting gas is assumed to be accelerated from rest at the 
inner edge of the disk, and thus is flowing at roughly the bal- 
listic free-fall speed at the stellar surface, 



Vff : 



2GAL 



^ R* 



1/2 



(8) 



This slightly underestimates the mean velocity, since in reality 
the streams come from radii between r m and r out . Using this 
expression allows the ram pressure at the stellar surface to be 
computed, with 



(9) 



(e.g., Hartmann et al. 1997; Calvet & Gullbring 1998). The 
accretion is assumed to be stopped at the point where P mm is 
balanced by the stellar atmosphere's gas pressure. A repre- 
sentative T Tauri model atmosphere (with gray opacity and 
local thermodynamic equilibrium) was used to determine the 
height dependence of density and gas pressure. The atmo- 
spheric density p s h was thus defined as that of the highest 
point in the atmosphere that remains undisturbed by either 
the shock or its (denser) post-shock cooling zone. The age 
dependence of p s h was computed from the condition of ram 
pressure balance. Examining these resulting values, though, 
yielded a useful approximation for this quantity. For most of 
the evolutionary ages considered for the solar-mass T Tauri 
star, the accretion is stopped a few scale heights above the 
photosphere, at which the temperature has reached its mini- 
mum radiative equilibrium value of — 0.87^. This value can 
be used in the definition of gas pressure to obtain a satisfac- 
tory estimate for p s h; i.e., one can solve for 



Ash 



5m H P ri 



(10) 



4k B T eS 

using equation (0 above and obtain a value within about 10% 
to 20% of the value interpolated from the full stellar atmo- 
sphere model. All of these values are only order-of-magnitude 
estimates, of course, since neither the post-shock heating nor 
the thermal equilibrium inside the accretion columns has been 
taken into account consistently. 

In addition to the density at the shock, the models below 
also require specifying the density at the visible stellar sur- 
face. A representative photospheric mass density was com- 
puted from the criterion that the Rosseland mean optical depth 
should have a value of approximately one: 



Tr f» K R p*H* = 1 



(11) 



where is the photospheric scale height as defined above 
and kr is the Rosseland mean opacity (in cm 2 g" 1 ) interpo- 
lated as a function of temperature and pressure from the table 
of Kurucz (1992). The resulting ratio p s h/ p* varies strongly 
over the T Tauri evolution, from a value of order unity (for 
the youngest stars with the strongest accretion) down to 10~ 3 
at the end of the accretion phase. 

Figure 4 shows several of the accretion-stream quantities 
defined above. Figure 4a shows the how the fractional area 
of both the accretion streams (5) and the polar caps (<5 po i) 
decrease similarly with increasing age. Calvet & Gullbring 
(1998) gave observationally determined values of 5 for a num- 
ber of T Tauri stars, but there was no overall trend with age; 
S was distributed (apparently) randomly between values of 
about 0.001 and 0. 1 for ages between log I0 1 = 4 and 7. Dupree 
et al. (2005) estimated <5 po i to be about 0.3, which agrees with 
the youngest T Tauri stars modeled here. Figure 4b plots the 
age dependence of both the photospheric density p* and the 
density at the accretion shock p s h (computed both numerically 
and using eq. IfTUI ). 

3.3. Properties of Inhomogeneous Accretion ("Clumps") 
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FIG. 4. — Accretion stream properties as a function of age. (a) Filling 
factor of accretion columns that thread the disk (solid line), filling factor of 
polar caps containing open field lines (dashed line), and the reciprocal of the 
number of inhomogeneous accretion tubes in one hemisphere (dotted line), 
(b) Mass densities at the stellar photosphere (solid li ne) a nd at the accretion 
shock (exact: dashed line, approximation from eq. 1101 : dotted line), and 
mean pre-shock density in the accretion stream (dot-dashed line). 

Observational evidence for intermittency and time variabil- 
ity in the accretion streams was discussed in § 2. Here, the 
flow along magnetospheric accretion columns is modeled as 
an idealized "string of pearls," where infalling dense clumps 
(i.e., blobs or clouds) are surrounded by ambient gas of much 
lower density. The clumps are assumed to be roughly spher- 
ical in shape with roughly the same extent in latitude as the 
magnetic flux tubes connecting the star and accretion disk. 
Thus, once they reach the stellar surface, the clumps have a 
radius given by 

^*(#in-f?out) 

r c = r (12) 

where the angles 6 m and 9 oul must be expressed in radians. 
Figure 2b shows how r c varies with age. 

The dense clumps are assumed to impact the star with a ve- 
locity v c equivalent to the free-fall speed vg defined in equa- 
tion (0. Scheurwater & Kuijpers (1988) defined a fiducial 
shock-crossing time for clumps impacting a stellar surface, 
with t c = \.5r c /v c . This time scale is of the same order of mag- 
nitude as the time it would take the blob to pass through the 
stellar surface if it was not stopped. It is useful to assume that 



the clumps are spread out along the flux tube with a constant 
periodicity; i.e., that along a given flux tube, clumps impact 
the star at a time interval At = Q c , where the dimensionless 
intermittency factor £ must be larger than 1 . 

When considering the mass accreted by infalling clumps, 
we follow Scheurwater & Kuijpers (1988) and define the den- 
sity interior to each clump as p c and the (lower) ambient inter- 
clump density as po. For the purpose of mass flux conserva- 
tion along the magnetic flux tubes, these densities are defined 
at the stellar surface. Despite the fact that the ratio p c /po ap- 
pears in several resulting quantities below, it always is can- 
celed out by other quantities that depend separately on p c and 
po such that the numerical value of their ratio never needs to 
be specified directly. A more relevant quantity for accretion 
is the mean density (p) in the flux tube, which is by definition 
larger than po and smaller than p c . The clump "overdensity ra- 
tio" (p c /(p)) can be computed by comparing the volume sub- 
tended by one clump with the volume traversed by a clump 
over time Af. In other words, if one assumes that p c ^S> po, 
one can find the mean density (p) by spreading out the gas in 
each clump to fill its own portion of the flux tube. Thus, 



Pc 

(p) 



r.At', 



47tr3/3 



(13) 



The overdensity ratio (or, equivalently, Q is a free parameter 
of this system, and a value of p c j (p) = 3 was chosen more 
or less arbitrarily. This sets £ = 24/9. Below, the resulting 
Alfven wave amplitude (caused by the periodic clump im- 
pacts) is found to depend on an overall factor of C}' 2 ; i.e., it 
is relatively insensitive to the chosen value of the overdensity 
ratio. 

In order to relate the inhomogeneous properties along a 
given flux tube to the total mass accretion rate M acc , the total 
number of flux tubes impacting the star must be calculated. 
The quantity is defined as the number of flux tubes in ei- 
ther the northern or southern hemisphere that contain accret- 
ing clumps. This definition is specific to the assumption of 
an aligned dipole magnetic field, and it is convenient because 
the summed effect of infalling clumps measured at the north 
[south] pole is assumed to depend only on the flux tubes in 
the northern [southern] hemisphere. The total number of flux 
tubes impacting the star is thus assumed to be 2N. One can 
compute by comparing two different ways of expressing the 
mass accretion rate. First, we know that in a time-averaged 
sense, mass is accreted with a local flux lp)v c over a subset of 
the stellar surface 5 given by equation (O. Thus, 



M acc = 4tt6RI( P )v c 



(14) 



Alternately, the mass in the 2N flux tubes can be summed up 
by knowing that each clump deposits a mass m c = 4iTr^p c /3, 
with 

2Nm c 16ttN , 



At 



9C 



■ PcV c r c 



(15) 



Equations < TT~4T > and (15[ must give the same total accretion 
rate, so equating them gives a useful expression for 



N = 2.) | y- 



(16) 



where equation ([T"3l was also used. This quantity is used be- 
low to compute the total effect of waves at the poles from the 
individual impact events. 
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Figure 4a shows N~ y (rather than N itself, to keep it on the 
same scale as the other plotted quantities), and it is interest- 
ing that N remains reasonably constant around 100-150 over 
most of the classical T Tauri phase of evolution. This mean 
value (for any instant of time) is large, but it is not so large 
that any fluctuations around the mean would be unresolvable 
due to averaging over the star. If the distribution of flux tubes 
is assumed to follow some kind of Poisson-like statistics, a 
standard deviation of order N 1 / 2 would be expected. In other 
words, for N « 100 there may always be something like a 10% 
level of fluctuations in the magnetospheric accretion rate. 

The equations above also allow (p) and p c to be computed 
from the known accretion rate M acc . Figure 4b shows (p) to 
typically be several orders of magnitude smaller than both p^ 
and p*. This large difference arises because the ram pressure 
inside the accretion stream is much larger than the gas pres- 
sure in the stream. The ratio (p) / p^ at the stellar surface is 
given very roughly by (c s /v c ) 2 , where c s is the sound speed 
corresponding to T e s, and c s <C v c . For the youngest T Tauri 
stars modeled, however, c s w v c and thus (p) « p S h. 

3.4. Properties of Impact-Generated Waves 

This section utilizes the results of Scheurwater & Kui- 
jpers (1988), who computed the flux of magnetohydrody- 
namic (MHD) waves that arise from the impact of a dense 
clump onto a stellar surface. Scheurwater & Kuijpers (1988) 
derived the total energy released in both Alfven and fast-mode 
MHD waves from such an impact in the "far-field" limit (i.e., 
at horizontal distances large compared to the clump size r c ). 
Slow-mode MHD waves were not considered because of the 
use of the cold-plasma approximation. The models below uti- 
lize only their Alfven wave result, since the efficiency of fast- 
mode wave generation was found to be significantly lower 
than for Alfven waves. Also, due to their compressibility, 
the magnetosonic (fast and slow mode) MHD waves are ex- 
pected to dissipate more rapidly than Alfven waves in stellar 
atmospheres (e.g., Kuperus et al. 1981; Narain & Ulmschnei- 
der 1990, 1996; Whang 1997). Thus, even if fast-mode and 
Alfven waves were generated in equal amounts at the impact 
site, the fast-mode waves may not survive their journey to 
the polar wind-generation regions of the star without being 
strongly damped. 

For simplicity, we assume that waves propagate away from 
the impact site with an overall energy given by the Scheurwa- 
ter & Kuijpers (1988) Alfven wave result, and that the waves 
undergo negligible dissipation. The wave energy released in 
one impact event is given by 



0.06615 



Pc 
IX) 



(17) 



where the numerical factor in front was given approximately 
as 0.07 in equation (57) of Scheurwater & Kuijpers (1988), 
but has been calculated more precisely from their Bessel- 
function integral. In this context, the Alfven speed V A is de- 
fined as that of the ambient medium, with 



Va 



Ba 



\J4npii 



(18) 



and Bq w 1000 G is the ambient magnetic field strength at the 
stellar surface, where the accretion streams connect with the 
star. Scheurwater & Kuijpers (1988) assumed that v c < Va 
(i.e., that the impacting clumps do not strongly distort the 
background magnetic field). This tends to lead to a very low 



efficiency of wave generation, such that equation (TTTb may be 
considered a "conservative" lower limit to the available en- 
ergy of fluctuations. For a given accretion column, the wave 
power that is emitted continuously by a stream of clumps is 
given by Ea/A?. 

In order to compute the wave energy density at other points 
on the stellar surface, we make the assumption that the waves 
propagate out horizontally from the impact point. It is impor- 
tant to note that ideal Alfven waves do not propagate any en- 
ergy perpendicularly to the background magnetic field. There 
are several reasons, however, why this does not disqualify the 
adopted treatment of wave energy propagation over the hori- 
zontal stellar surface. First, the true evolution of the impact 
pulse is probably nonlinear. Much like the solar Moreton and 
EIT waves discussed in § 2, nonlinear effects such as mode 
coupling, shock steepening, and soliton-like coherence are 
likely to be acting to help convey the total "fluctuation yield" 
of the impact across the stellar surface. Second, for any real 
star, the background magnetic field is never completely radial, 
and it will always have some component horizontal to the sur- 
face (along which even linear Alfven waves can propagate). 
Thus, the dominant end-result of multiple cloud impacts is 
assumed here to be some kind of transverse field-line pertur- 
bations that are treated for simplicity with the energetics of 
Alfven waves. 

Considering waves that spread out in circular ripples from 
the impact point, the goal is to compute the horizontal wave 
flux (power per unit area) at a distance x away from the central 
point. For this purpose, the stellar surface can be treated as 
a flat plane. The wave power is assumed to be emitted into 
an expanding cylinder with an approximate height of 2r c (the 
clump diameter) and a horizontal radius x. The horizontal flux 
Fa of waves into the vertical wall of the cylinder is given by 
dividing the power by the area of the surrounding wall, with 



£a/At 
4nxr r 



0.0147 



2 3 



x Po( \Va 



(19) 



Note that the flux decreases linearly with increasing x, as is 
expected for a cylindrical geometry. 

The total wave flux from the effect of multiple clump im- 
pacts is computed at the north or south pole of the star. The 
accretion stream impact points are assumed to be distributed 
circularly in a ring around the pole. Thus, using the geomet- 
ric quantities derived earlier, this implies x = R*{9i n + 9 out )/2. 
Also, because each accretion column is assumed to be at the 
same horizontal distance from the target point at the pole, the 
total wave flux is given straightforwardly by N times the flux 
due to one stream of clumps. The total accretion-driven wave 
flux arriving at either the north or south pole is thus given by 
Fuot = NFa- It is assumed that the waves do not damp appre- 
ciably between where they are generated (at the bases of the 
accretion streams) and their destination (at the pole). 

The standard definition of the flux of Alfven waves, in a 
medium where the bulk flow speed is negligible in compari- 
son to the Alfven speed, is 



Fam = PoVjVa 



(20) 



This can be compared with the total wave flux derived above 
to obtain the perpendicular Alfven wave velocity amplitude 
v_i_. In units of the clump infall velocity v c , the wave ampli- 
tude is thus given by 

Zi = 0.1715 P-l(^-\' 2 (^] 2 . (21) 
Vc Po V & J 



Va 



8 



S. R. CRANMER 



Note that there is no actual dependence on po, since the ex- 
plicit factor of po in the denominator is canceled by the den- 
sity dependence of V^ 2 . 

The wave amplitude derived in equation (I2TI 1 is the value at 
the shock impact height, which formally can be above or be- 
low the photosphere. The stellar wind models below, though, 
require the Alfven wave amplitude to be specified exactly at 
the photosphere. If we assume that the wave energy density 
Ua = P v2 ± is conserved between the shock height and the pho- 
tosphere, then the densities at those heights can be used to 
scale one to the other, with 



10" 



V_L 



Ash 
P* 



1/2 



(22) 



where vj_* is the photospheric wave amplitude, and the densi- 
ties p* and p s h were defined above in § 3.2. 

Although the overall energy budget of accretion-driven 
waves is treated under the assumption that the waves are 
Alfvenic, it is likely that the strongly nonlinear stream impacts 
also give rise to compressible waves of some kind. As men- 
tioned above, the analysis of Scheurwater & Kuijpers (1988) 
did not take account of slow-mode MHD waves that, for par- 
allel propagation and a strong background field, are identical 
to hydrodynamic acoustic waves. There is evidence, however, 
that another highly nonlinear MHD phenomenon — turbulent 
subsurface convection — gives rise to both longitudinal (com- 
pressible) and transverse (incompressible) MHD waves with 
roughly comparable energy densities (e.g., Musielak & Ulm- 
schneider 2002). Thus, the models below are given a photo- 
spheric source of accretion-driven acoustic waves that are in 
energy equipartition with the accretion-driven Alfven waves; 
i.e., U s = Ua- The upward flux of acoustic waves is thus given 
by F s = c s U s , where c s is the sound speed appropriate for T e ff. 

For both the acoustic and Alfven waves at the photo- 
sphere, the accretion-driven component is added to the intrin- 
sic (convection-driven) component. A key assumption of this 
paper is that the convection-driven component is held fixed, as 
a function of age, at the values used by Cranmer et al. (2007) 
for the present-day Sun. This results in minimum values for 
the Alfven wave amplitude (0.255 km s" 1 ) and the acoustic 
wave flux (10 8 erg cm" 2 s -1 ) below which the models never 
go. There are hints that rapidly rotating young stars may un- 
dergo more intense subsurface convection than the evolved 
slowly rotating Sun (Kapyla et al. 2007; Brown et al. 2007; 
Ballot et al. 2007), but the implications of these additional 
variations with age are left for future work (see also § 6). 

Figure 5 displays various velocity quantities used in the 
accretion-driven wave scenario. The clump free-fall speed 
v c = Vff always remains smaller than the ambient Alfven speed 
Va, which was an assumption that Scheurwater & Kuijpers 
(1988) had to make in order for the background magnetic 
field to remain relatively undisturbed by the clumps. The 
accretion-driven Alfven wave amplitude at the shock is larger 
than that measured at the photosphere (see eq. El ), and at 
ages later than about 60 Myr the accretion-driven waves at 
the photosphere grow weak enough to be overwhelmed by the 
convection-driven waves. For a limited range of younger ages 
(2 x 10 4 < t < 5 x 10 5 yr) the Alfven wave motions are super- 
sonic in the photosphere, with a peak value of v_u = 18 km 
s -1 at t = 5 x 10 4 yr. The sharp decrease in wave amplitude at 
the youngest ages is due to the inner edge of the accretion disk 
coming closer to the stellar surface. In that limit, the ballistic 
infall speed Vft becomes small and the latitudinal distance x 
traversed by the waves becomes large, thus leaving negligible 
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FIG. 5. — Velocities related to accretion streams as a function of age: free- 
fall clump speed (thin solid line), ambient Alfven speed (dot-dashed line), 
photospheric sound speed (dotted line). Plotted Alfven wave amplitudes are 
those due to accretion-streams and measured at the shock (dashed line) and 
those due to both accretion and convection, measured at the photosphere 
(thick solid line). 

energy in the waves once they reach the pole. 

Finally, it is worthwhile to sum up how the above values 
depend on the relatively large number of assumptions made 
about the accretion streams. First, the use of an aligned dipole 
field for the accretion streams, the Konigl (1991) expression 
for r m , and the "string of pearls" model of clumps along the 
flux tubes are all somewhat simplistic and should be replaced 
by more realistic conditions in future work. Second, there 
are three primary parameters that had to be specified in or- 
der to determine numerical values for the accretion properties. 
These are: (1) the relative size of the outer disk radius with re- 
spect to the inner disk radius; i.e., e = 0.1, (2) the clump over- 
density ratio p c j (p) = 3, and (3) the magnetic field strength at 
the base of the accretion streams Bo = 1000 G. Third, proba- 
bly the most idealistic assumption in the modeled evolution- 
ary sequence is the use of a single monotonic relation for M acc 
versus t (e.g., eq. Q). § 6 describes how these assumptions 
can be relaxed in subsequent modeling efforts. 

4. IMPLEMENTATION IN STELLAR WIND MODELS 

The steady-state outflow models presented here are numer- 
ical solutions to one-fluid conservation equations for mass, 
momentum, and energy along a polar flux tube. For the spe- 
cific case of the solar wind, Cranmer et al. (2007) presented 
these equations and described their self-consistent numeri- 
cal solution using a computer code called ZEPHYR. The T 
Tauri wind models are calculated with an updated version of 
ZEPHYR, with specific differences from the solar case de- 
scribed below. 

4.1. Conservation Equations and Input Physics 

The equation of mass conservation along a magnetic flux 
tube is 

1 d 

- — ( P uA) = (23) 
A or 

where u and p are the wind speed and mass density specified 
as functions of radial distance r, and A is the cross-sectional 
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area of the polar flux tube. Magnetic flux conservation de- 
mands that the product BoA is constant along the flux tube, 
where Brj(0 is the field strength that is specified explicitly. 
The equation of momentum conservation is 



du 1 dP 
dr p dr 



GM* 



+ D 



(24) 



where P is the gas pressure and D is the bulk acceleration 
on the plasma due to wave pressure; i.e., the nondissipative 
net ponderomotive force due to wave propagation through the 
inhomogeneous medium (Bretherton & Garrett 1968; Belcher 
1971; Jacques 1977). A complete expression for D in the 
presence of damped acoustic and Alfven waves was given by 
Cranmer et al. (2007). The simpler limit of wave pressure due 
to dissipationless Alfven waves is discussed in more detail in 
§§ 5.2-5.3. 

For the pure hydrogen plasma assumed here, the gas pres- 
sure is given by P = (1 +x)pksT /run, where x is the hydro- 
gen ionization fraction. Note that although the pressure is 
calculated for a hydrogen gas, the radiative cooling rate Q m & 
used in the energy equation is dominated by metals. Cranmer 
et al. (2007) tested the ZEPHYR code with two separate as- 
sumptions for the ionization balance. First, a self-consistent, 
but computationally intensive solution was used, which im- 
plemented a three-level hydrogen atom. In that model, the 
n = 1 and n = 2 levels were assumed to remain in relative local 
thermodynamic equilibrium (LTE) and the full rate equation 
between n = 2 and the continuum was solved iteratively. Sec- 
ond, a simpler tabulated function of x as a function of tem- 
perature T was taken from a semi-empirical non-LTE model 
of the solar photosphere, chromosphere, and transition region 
(e.g., Avrett & Loeser 2008). For both the solar and T Tauri 
star applications, the results for the two cases were extremely 
similar, and thus the simpler tabulated function was used in 
the models described below. 

For solar-type winds, the key equation for both heating the 
corona and setting the mass loss rate is the conservation of 
internal energy, 

BE fE+P\ d 

U—+ [ — 1 TriuA) = grad+econd + fiA + eS (25) 



dr 



dr 



where E is the internal energy density and the terms on the 
right-hand side are volumetric heating/cooling rates due to 
radiation, conduction, Alfven wave damping, and acoustic 
(sound) wave damping. The terms on the left-hand side that 
depend on u are responsible for enthalpy transport and adi- 
abatic cooling. In a partially ionized plasma, the definition 
of E is convention-dependent; we use the same one as Ulm- 
schneider & Muchmore (1986) and Mullan & Cheng (1993), 
which is E = (3P /2) + xpln/ mjj, where In = 13.6 eV. The net 
heating/cooling from conduction (g con d) is given by a gradual 
transition between classical Spitzer-Harm conductivity (when 
the electron collision rate is fast compared to the wind ex- 
pansion rate) and Hollweg's (1974, 1976) prescription for 
free-streaming heat flux in the collisionless heliosphere (when 
electron collisions are negligible). All of these terms were de- 
scribed in more detail by Cranmer et al. (2007). 

The volumetric radiative heating/cooling rate Q m & has been 
modified slightly from the earlier solar models. The solar rate 
is first computed as in Cranmer et al. (2007), but then it is 

4 Although this is formally inconsistent, the resulting properties of the 
plasma are not expected to be far from values computed with a more accurate 
equation of state. 



multiplied by a correction factor similar to that suggested by 
Hartmann & Macgregor (1980) for massive, optically thick 
chromospheres of late-type stars. The correction factor / is 
assumed to be proportional to the escape probability P esc for 
photons in the core of the Mg II A2803 (h) and A2796 (k) lines. 
A simple expression that bridges the optically thin and thick 
limits, for a line with Voigt wings, is 

/ = 2P esc w —Ljf (26) 

1 +T ' 

(e.g., Mihalas 1978), which may err on the side of overesti- 
mating the escape probability and thus would give a conserva- 
tive undercorrection to Q md . Hartmann & Macgregor (1980) 
assumed that the optical depth of low-ionization metals would 
scale as g~^ 2 , where g is the stellar surface gravity. Here, we 
compute the optical depth in the core of the h and k lines (Thk) 
more exactly by integrating over the radial grid of density and 
temperature in the wind model, 



= / Xhk v7 



dr . 



(27) 



where we include the spherical correction factor suggested by 
Lucy (1971, 1976) for extended atmospheres. The line-center 
extinction coefficient is given approximately by 



Xhk = 0.0153 



/ »Mg 



VlO 10 



cm J / 



10 4 K 



-1/2 



cm 



(28) 



and the number density «Mgn of ground-state ions is given by 
the product of the hydrogen number density (p/m H ), the Mg 
abundance with respect to hydrogen (3.4 x 10" 5 ; Grevesse et 
al. 2007), and the ionization fraction of Mg + with respect to 
all Mg species. For temperatures in excess of 10 4 K, the lat- 
ter is given by coronal ionization equilibrium (Mazzotta et al. 
1998); for temperatures below 10 4 K, the coronal equilibrium 
curve is matched onto a relation determined from LTE Saha 
ionization balance. 

Following Hartmann & Macgregor (1980), the opacity cor- 
rection factor / is ramped up gradually as a function of tem- 
perature. When T < 7160 K, Q m ^ is multiplied by the full 
factor /. For 7160 < T < 8440 K, the multiplier is varied by 
taking y( 844 °- 7 ')/ 1280 5 which ramps down to unity by the time 
the temperature rises to 8440 K. For temperatures in excess of 
8440 K, the correction factor is not used. 

The evolution and damping of Alfven and acoustic waves — 
which affects both D in the momentum equation and Q A and 
Qs in the energy equation — is described in much more detail 
by Cranmer et al. (2007). Given a specified frequency spec- 
trum of acoustic and Alfven wave power in the photosphere, 
equations of wave action conservation are solved to determine 
the radial evolution of the wave amplitudes. Acoustic waves 
are damped by both heat conduction and entropy gain at shock 
discontinuities. Alfven waves are damped by MHD turbu- 
lence, for which we only specify the net transport of energy 
from large to small eddies, assuming the cascade must termi- 
nate in an irreversible conversion of wave energy to heat. 

Coupled with the wave action equations are also non-WKB 
transport equations to determine the degree of linear reflection 
of the Alfven waves (e.g., Heinemann & Olbert 1980). This 
is required because the turbulent dissipation rate depends on 
differences in energy density between upward and downward 
traveling waves (see also Matthaeus et al. 1999; Dmitruk et al. 
2001, 2002; Cranmer & van Ballegooijen 2005). The result- 
ing values of the Elsasser amplitudes Z±, which denote the 



10 



S. R. CRANMER 



energy contained in upward (Z_) and downward (Z+) propa- 
gating waves, were then used to constrain the energy flux in 
the cascade. The Alfven wave heating rate (erg s" 1 cm" 3 ) is 
given by a phenomenological relation that has evolved from 
analytic studies and numerical simulations; i.e., 

e »='(TTwd-«r^ ,29) 

(see also Hossain et al. 1995; Zhou & Matthaeus 1990; 
Oughton et al. 2006). The transverse length scale £± is an 
effective perpendicular correlation length of the turbulence, 
and Cranmer et al. (2007) used a standard assumption that £± 
scales with the cross-sectional width of the flux tube (Holl- 
weg 1986). The term in parentheses in equation (1291 is an ef- 
ficiency factor that accounts for situations in which the turbu- 
lent cascade does not have time to develop before the waves or 
the wind carry away the energy (Dmitruk & Matthaeus 2003). 
The cascade is quenched when the nonlinear eddy time scale 
f e ddy becomes much longer than the macroscopic wave reflec- 
tion time scale f re f . 

4.2. Model Inputs and Numerical Procedures 

The ZEPHYR code was designed to utilize as few free pa- 
rameters as possible. For example, the coronal heating rate 
and the spatial length scales of wave dissipation are computed 
internally from straightforward physical principles and are 
not input as adjustable parameters. However, there are quanti- 
ties that do need to be specified prior to solving the equations 
given in § 4. 1 . 

An important input parameter to ZEPHYR is the radial de- 
pendence of the background magnetic field Bo(r). Although 
many stellar field-strength measurements have been made, 
relatively little is known about how rapidly Bo decreases with 
increasing height above the photosphere, or how fragmented 
the flux tubes become on the granulation scale. Some in- 
formation about this kind of structure is contained in "fill- 
ing factors" that can be determined observationally (e.g., Saar 
2001). However, these measurements may be biased toward 
the bright closed-field active regions and not the footpoints of 
stellar wind streams. Thus, as in other areas of this study, the 
present-day solar case was adopted as a baseline on which to 
vary the evolving properties of the 1 M model star. 

The multi-scale polar coronal hole field of Cranmer & van 
Ballegooijen (2005) was used as the starting point for the 
younger T Tauri models. The lower regions of that model 
were derived for a magnetostatic balance between the strong- 
field (low density) solar wind flux tube and the weak-field 
(high density) surrounding atmosphere. The radial depen- 
dence of magnetic pressure Bq/8tt thus scales with the gas 
pressure P. In the T Tauri models, then, the lower part of 
the magnetic field was stretched in radius in proportion with 
the pressure scale height Figure 6 shows Bo(r) for several 
ages, along with an idealized dipole that does not take account 
of the lateral expansion of the flux tube close to the surface. 

In addition to the global magnetic field strength, the 
ZEPHYR models also require three key wave-driving param- 
eters to be specified at the lower boundary: 

1 . The photospheric acoustic flux F s mainly affects the 
heating at chromospheric temperatures (T ~ 10 4 K). 
The solar value of 10 8 erg cm" 2 s" 1 was summed with 
the accretion-driven acoustic flux as discussed in § 3.4. 
The power spectrum of acoustic waves was adopted 
from the solar spectrum given in Figure 3 of Cranmer 
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FIG. 6. — Background magnetic field Bo over the stellar poles as a function 
of the height above the photosphere (measured in stellar radii). Present-day 
solar field strength (thick solid line) is compared with that of ages: log to t = 
5.5 (dotted line), 6.0 (dot-dashed line), and 6.5 (dashed line). Also shown is 
an ideal dipolar field having the same strength at the photosphere as the other 
cases (thin solid line). 

et al. (2007), and the frequency scale was shifted up or 
down with the photospheric value of the acoustic cutoff 
frequency, oj. dc = c s j2H % , which evolves over time. 

2. The photospheric Alfven wave amplitude v^* is speci- 
fied instead of the upward energy flux Fa because the 
latter depends on the cancellation between upward and 
downward propagating waves that is determined as a 
part of the self-consistent solution. As discussed in 
§ 3.4, the solar value of 0.255 km s" 1 was supplemented 
by the accretion-driven wave component. The shape of 
the Alfven wave frequency spectrum was kept fixed at 
the solar model because it is unclear how it should scale 
with varying stellar properties. 

3. The photospheric Alfven wave correlation length £j_* 
sets the scale of the turbulent heating rate Qa (eq. l29l ). 
Once this parameter is set, the value of £± at larger 
distances is determined by the assumed proportional- 
ity with A 1 / 2 . The solar value of 75 km (Cranmer et 
al. 2007) was evolved in proportion with changes in the 
photospheric scale height The justification for this 
is that the horizontal scale of convective granulation is 
believed to be set by the scale height (e.g., Robinson 
et al. 2004), and the turbulent mixing scale is probably 
related closely to the properties of the granulation. 

The numerical relaxation method used by ZEPHYR was 
discussed by Cranmer et al. (2007). In the absence of ex- 
plicit specification here, the new T Tauri wind models use the 
same parameters as given in that paper. However, the new 
models use a slightly stronger form for the code's iterative 
undercorrection than did the original solar models. From one 
iteration to the next, ZEPHYR replaces old values with with 
a fractional step toward the newly computed values, rather 
than using the new values themselves. This technique was 
motivated by globally convergent backtracking methods for 
finding roots of nonlinear equations (e.g., Dennis & Schnabel 
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1983). The solar models used a constant minimum under- 
correction exponent eo = 0.17, as defined in equation (65) of 
Cranmer et al. (2007). The T Tauri models started with this 
value, but gradually decreased it over time by multiplying eo 
by 0.97 after each iteration. The value of eo was not allowed 
to be smaller than an absolute minimum of 0.001. This rep- 
resented an additional kind of "annealing" that helped the pa- 
rameters reach their time-steady values more rapidly and ro- 
bustly. 

5. RESULTS 
5.1. Standard Age Sequence 

A series of 24 ZEPHYR models was created with ages 
ranging between log 10 f = 5.65 and 9.66. These models all 
converged to steady-state mass-conserving wind outflow so- 
lutions within 250 numerical iterations. The relative errors in 
the energy equation, averaged over the radial grid, were of or- 
der 1 % for the final converged models. For ages younger than 
log 1() ? = 5.65 (i.e., 0.45 Myr), it was found that time-steady 
solutions to equations d23ll-(l2"5ll do not exist, and the best we 
can do is to estimate the mass flux that is driven up through the 
Parker critical point (see §§ 5.2-5.3 below). Figure 7 presents 
a summary of various scalar properties for the wind models as 
a function of age. 

The mass loss rate M w ; n( j, shown in Figure la, was calcu- 
lated by multiplying the mass flux pu at the largest grid radius 
(1200/?* for all models) by the full spherical area 4-nr 2 at that 
radius. This is a slight overestimate, since the polar flux tubes 
only cover a finite fraction of the stellar surface {S po i < 1). 
However, it is unknown whether these regions expand out- 
ward or inward (i.e., to larger or smaller solid angle coverage) 
with increasing distance from the star. The actual large-scale 
geometry is likely to depend on how the pressure (gas and 
magnetic) in the field lines that thread the accretion disk is 
able to confine or collimate the polar flux tubes. 

As the age is decreased, M w i n d for the time-steady mod- 
els increases by four orders of magnitude, from the present- 
day solar value of about 2 x 10" 14 M Q yr" 1 up to 4 x 10" 10 
Mq yr" 1 at the youngest modeled age of 0.45 Myr. Note that 
these values exceed the mass loss rates predicted by the empir- 
ical scaling relation of Reimers (1975, 1977) at all ages (i.e., 
using M w ; n[ j oc L^Rx/M* and normalizing to the present-day 
mass loss rate), but for t > 100 Myr there is rough agreement. 
The stellar wind velocity at the largest radius is denoted as 
an "asymptotic" or terminal speed and is shown in Fig- 
ure lb. The wind speed remains roughly constant for most of 
the later phase of the evolution (with « V e scX but it drops 
precipitously in the youngest models. 

The dominant physical processes that drive the evolving 
stellar winds are revealed when examining the maximum tem- 
peratures in the models, as shown in Figure 7c. The older 
(t > 1 .75 Myr) models with high wind speeds and solar-like 
mass loss rates have hot coronae, with peak temperatures be- 
tween 10 6 and 2 x 10 6 K. The younger models undergo a 
rapid drop in temperature, ultimately leading to "extended 
chromospheres" with peak temperatures around 10 4 K. This 
transition occurs because of a well-known thermal instabil- 
ity in the radiative cooling rate 2 ra d (e.g., Parker 1953; Field 
1965; Rosner et al. 1978; Suzuki 2007). At temperatures be- 
low about 10 5 K, the cooling rate decreases with decreasing 
T, and small temperature perturbations are easily stabilized. 
However, above 10 5 K, Qrad decreases as T increases, which 
gives any small increase in temperature an unstable runaway 
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FIG. 7. — Age-dependent model properties, (a) Mass loss rates of time- 
steady ZEPHYR models (solid line with symbols) compared with the Reimers 
(1975, 1977) relation (dashed line), analytic/cold models of § 5.3 (dotted 
line), and adopted accretion rate (dot-dashed line), (b) Terminal wind speed 
for ZEPHYR models (thick solid line with symbols), surface escape speed 
(dashed line), photospheric sound speed (dot-dashed line), and wind speed 
at the wave-modified critical point for ZEPHYR models (thin solid line) and 
analytic/cold models (dotted line), (c) Peak temperatures of ZEPHYR mod- 
els (solid line with symbols), temperatures at critical point for analytic/cold 
models (dotted line), and stellar r c ff (dot-dashed line). 

toward larger values. This is the same instability that helps 
lead to the sharp transition region between the 10 4 K solar 
chromosphere and the 10 6 K corona (see also Hammer 1982; 
Withbroe 1988; Owocki 2004). 

The reason that the young and old models in Figure 7 end 
up on opposite sides of the thermal instability is that the to- 
tal rates of energy input (i.e., the wave heating rates Q A and 
<2s) vary strongly as a function of age. The older models have 
lower wave amplitudes and thus weaker heating rates, which 
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leads to relatively low values of M W md- The correspondingly 
low atmospheric densities in these models give rise to weak 
radiative cooling (because oc p 2 ) that cannot suppress the 
coronal heating. The coronal winds are driven by comparable 
contributions from gas pressure and wave pressure. On the 
other hand, the younger models have larger wave amplitudes, 
more energy input (as well as more wave pressure, which ex- 
pands the density scale height), and thus more massive winds 
with stronger radiative cooling. These chromospheric winds 
are driven mainly by wave pressure. Related explanations for 
cool winds from low-gravity stars have been discussed by, 
e.g., Antiochos et al. (1986), Rosner et al. (1991, 1995), Will- 
son (2000), Killie (2002), and Suzuki (2007). 

The processes that determine the quantitative value of M wm d 
are different on either side of the thermal bifurcation. The 
mass flux for the older, solar-like models can be explained 
well by the energy balance between heat conduction, radia- 
tive losses, and the upward enthalpy flux. The atmosphere 
finds the height of the transition region that best matches these 
sources and sinks of energy in a time-steady way, and the re- 
sulting gas pressure at this height sets the mass flux. Various 
analytic solutions for this balance have been given by Ham- 
mer (1982), Withbroe (1988), Leer et al. (1998), and Cranmer 
(2004). The younger, more massive wind models are found 
to approach the limit of "cold" wave-driven outflows, where 
the Alfven wave pressure replaces the gas pressure as the pri- 
mary means of canceling out the stellar gravity. In this limit, 
the analytic approach of Holzer et al. (1983), discussed fur- 
ther in § 5.3, has been shown to provide a relatively simple 
estimate for the mass loss rate. 

Figure 8 illustrates a selection of other parameters of the 
ZEPHYR models. The heights of the wave-modified criti- 
cal point (see eq. l32l below) and the Alfven point (where 
// = Va) are shown in units of the stellar radius. The ratio of 
the isothermal sound speed a to the wind speed at the critical 
point is also shown. This ratio is close to unity for the older, 
less massive models (indicating that gas pressure dominates 
the wind acceleration) and is less than 1 for the younger (more 
wave-driven) models. Note also that the older "coronal" mod- 
els have critical points in the sub-Alfvenic wind (u < Va) 
and the younger "extended chromosphere" models have larger 
critical radii that are in the super- Alfvenic (u > Va) part of the 
wind. Also plotted in Figure 8a is the so-called wind lumi- 
nosity L w i n d, which we estimate from the sum of the energy 
required to lift the wind out of the gravitational potential and 
the remaining kinetic energy of the flow, i.e., 

,2' 



^wind — ^^wind 



u 
+ — 

2 



(30) 



(e.g., Clayton 1966). This expression ignores thermal energy, 
ionization energy, magnetic energy, and waves, which are ex- 
pected to be small contributors at and above the critical point. 
Curves for the ratio L w i„d /£* are shown both for the wind at 
its largest height (i.e., using Uoa in eq. l30l ) and at the critical 
point (using « cr i t ). The latter is helpful to compare with ana- 
lytic estimates for younger ages described below — for which 
nothing above the critical point was computed. 

Figure %b shows the terms that dominate the energy con- 
servation equation (at the critical point) as a function of age. 
The areas plotted here were computed by normalizing the ab- 
solute values of the individual terms in equation ( f25l l by the 
maximum absolute value for each model and then "stacking" 
them so that together they fill the region between and 1. 
The older coronal models have a three-part balance between 
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FIG. 8. — Additional age-dependent model properties, (a) Wave-modified 
critical radii and Alfven radii {upper solid lines, labeled), ratio of sound speed 
to wind speed at critical radius {dot-dashed line), and ratio of wind luminos- 
ity to photon luminosity, measured at the largest modeled radius {lower solid 
line) and at the critical point {dashed line). Results for the analytic/cold mod- 
els are also shown {dotted lines), {h) Areas denote the relative contribution 
of terms to the energy conservation equation at the critical point; see labels. 
(2adv denotes the advection terms on the left-hand side of equation (25), other 
terms are denned in the text.) Also plotted is the exponent £ versus age {dot- 
dashed line). 



Alfven wave heating, heat conduction, and the upward ad- 
vection of enthalpy due to the terms on the left-hand side of 
equation ( 1251 ). For the younger models, radiative losses be- 
come important because of the higher densities at the critical 
point, and heat conduction disappears because the tempera- 
tures are so low. 

Figure 9 displays the radial dependence of wind speed, tem- 
perature, and density for three selected ages: (1) the present- 
day polar outflow of the Sun, (2) a younger, but still low-M w j n d 
coronal wind, and (3) an even younger model that has made 
the transition to a higher mass loss rate and an extended chro- 
mosphere. The plotted hydrogen number density tin is that 
of all hydrogen nuclei (both neutral and ionized). Note that 
the model with the lowest temperature would be expected to 
have a small density scale height (eq. [H) an d thus a rapid ra- 
dial decline in «h- However, Figure 9c shows the opposite to 
be the case: the coolest model has the largest effective density 
scale height of the three because of a much larger contribution 
from wave pressure. 
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FIG. 9. — Radial dependence of wind parameters for three selected 
ZEPHYR models: (a) wind outflow velocity, (b) temperature, and (c) to- 
tal hydrogen number density. In all panels, models are shown for ages 
log 10 f = 6.0 (dotted lines), 6.5 (dashed lines), and 9.66 (solid lines). Wave- 
modified critical radii (circles) and Alfven radii (triangles) are also shown. 

There are two additional points of comparison with other 
work that should be made in the light of the main results 
shown in Figures 7-9. 

1. There is a relatively flat age dependence of the pre- 
dicted mass loss rate for the post-T Tauri phase (i.e., 
zero-age main sequence [ZAMS] stars with t > 50 
Myr). This stands in contrast to the observationally in- 
ferred power-law decrease for these stages, which is ap- 
proximately M w i n d oc f~ 2 (e.g., Wood et al. 2002, 2005). 
Note, however, that the internal (convection-related) 
source of MHD wave energy at the photosphere was 
assumed in the models to be set at the present-day solar 
level and not to vary with age. As described in § 6, the 
inclusion of rotation-dependent convection may give 



rise to a stronger variation in the turbulent fluctuation 
amplitude with age, and thus also a more pronounced 
age dependence of the rates of coronal heating and mass 
loss. 

2. Suzuki (2007) modeled the extended atmospheres and 
winds of low-gravity cool stars, and found that the on- 
set of thermal instability gave rise to large amounts of 
time variability. These models often showed dense and 
cool shells that coexist with spontaneously produced 
hot and tenuous "bubbles." Suzuki (2007) noted that 
this dynamical instability began to occur once the es- 
cape velocity V esc (measured in the wind acceleration 
region) dropped to the point of being of the same order 
of magnitude as the sound speed corresponding to the 
thermal instability temperature of ~ 10 5 K. It is unclear, 
though, whether this variability is triggered only for rel- 
atively moderate rates of turbulent energy input, like the 
amplitudes derived by Suzuki (2007) from non-rotating 
convection models. The larger amplitudes used in the 
present models of accretion-driven waves may be suf- 
ficient to drive the winds much more robustly past the 
thermal instability and into time-steady extended chro- 
mospheres (see also Killie 2002). 

5.2. Disappearance of Time-Steady Solutions 

The ZEPHYR code could not find time-independent wind 
solutions for ages younger than about 0.45 Myr (i.e., for ac- 
cretion rates larger than about 7 x 10~ 8 M yr" 1 ). This was 
found not to be a numerical effect. Instead, it is an outcome 
of the requirement that time-steady winds must have a suffi- 
cient amount of outward acceleration (either due to gas pres- 
sure or some other external forcing, like waves) to drive ma- 
terial out of the star's gravitational potential well to a finite 
coasting speed at infinity. This was realized for polytropic 
gas-pressure-driven winds very early on (Parker 1963; Holzer 
& Axford 1970), such that if the temperature decreases too 
rapidly with increasing distance (i.e., with decreasing density) 
the wind would become "stalled" and not have a time-steady 
solution. 

To illustrate how this effect occurs for the youngest, wave- 
driven models, let us write and analyze an approximate equa- 
tion of momentum conservation. For simplicity, the wind 
temperature T is assumed to be constant, and the radii of in- 
terest are far enough from the star such that the flux-tube ex- 
pansion can be assumed to be spherical, with A oc r 2 . Also, 
for these models the acoustic wave pressure can be ignored 
(since compressive waves damp out rather low in the atmo- 
sphere), and the radial behavior of the Alfven wave amplitude 
can be modeled roughly in the dissipationless limit. Thus, one 
can follow Jacques (1977) and write the momentum equation 
as a modified critical point equation 

\ u ) dr r l r 

At the wave-modified critical point (r CI i t ), the wind speed u 
equals the critical speed u a n, which is defined as 

2 2 v 2 /1 + 3M A \ 

"crit = a (~L+Ma~ J ^ 

where the squared isothermal sound speed is a 2 = (1 + 
x)k B T /m H , and the bulk-flow Alfven Mach number M A = 
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u/Va- A more general version of equation OTI ) that also con- 
tains damping and acoustic wave pressure is given in, e.g., § 6 
of Cranmer et al. (2007). 

The above expressions show how the outward pressure, 
which balances gravity at the critical point, can be dominated 
either by a 2 (gas pressure) or by a term proportional to v\ 
(wave pressure). For gas pressure that can be described as 
a polytrope (i.e., a 2 cx p 7_1 ), the polytropic index 7 at the 
critical point must be smaller than 1.5 in order for there to 
be a time-steady acceleration from the critical point to infin- 
ity (Parker 1963; see also Velli 2001; Owocki 2004). Larger 
polytropic indices 7 > 1.5 imply that a 2 drops too rapidly with 
increasing distance to provide sufficient acceleration for a par- 
cel accelerated through the critical point to escape to infinity. 

For winds dominated by wave pressure, it is possible to use 
the equation of wave action conservation (eq. l36l ) to exam- 
ine the density dependence of the wave amplitude in a sim- 
ilar way as above (e.g., Jacques 1977; Heinemann & Olbert 
1980; Cranmer & van Ballegooijen 2005). The exponent £ in 
the scaling relation vj_ oc p* is known to be a slowly varying 
function of distance. Close to the star, where u <C Va, the ex- 
ponent £ w -0.25. In the vicinity of the Alfven point (u « Va), 
£ increases to zero and v_i_(r) has a local maximum. Far from 
the star, where u Va, the exponent £ grows to an asymp- 
totically large value of +0.25. The dimensionless quantity in 
parentheses in equation (l32l varies only between 1 and 3 over 
the full range of distances, and can be assumed to be roughly 
constant compared to the density dependence of vj_. 

Thus, by comparing the polytropic, gas-pressure dominated 
expression for M 2 . it to the wave-pressure dominated version, it 
becomes possible to write an effective polytropic index for the 
latter case as 7 e ff = 2£+ 1 . The unstable region of 7 > 1 .5 cor- 
responds to £ > 0.25, and this value is reached at the critical 
point only when r cl -i t is well into the super- Alfvenic part of the 
wind. Indeed, this corresponds to the youngest models shown 
in Figures 7-9. As M w ; n d increases (with decreasing age), the 
density in the wind increases, and this leads to a sharp decline 
in the value of Va at the critical point. (Just over the span of 
ages going from log 10 f = 6.25 to 6.0, Va at the critical point 
drops by a factor of 150.) Figure 8Z? shows £ versus age for 
the ZEPHYR models, where this exponent was first computed 
for all heights via 



d(\nv±)/dr 



(33) 



d(\np)/dr 

and the value shown is that interpolated to the location of 
the wave-modified critical point. It is clear that £ approaches 
the unstable limiting value of 0.25 just at the point where the 
time-steady solutions disappear. 

What happens to a stellar outflow when "too much" mass 
is driven up past the critical point to maintain a time-steady 
wind? An isolated parcel of gas with an outflow speed u = w cr i t 
at the critical point would be decelerated to stagnation at some 
height above the critical point, and it would want to fall back 
down towards the star. In reality, this parcel would collide 
with other parcels that are still accelerating, and a stochastic 
collection of shocked clumps is likely to result. Interactions 
between these parcels may result in an extra degree of col- 
lisional heating that could act as an extended source of gas 
pressure to help maintain a mean net outward flow. Situations 
similar to this have been suggested to occur in the outflows 
of both pulsating cool stars (e.g., Bowen 1988; Willson 2000; 
Struck et al. 2004) and luminous blue variables (Owocki & 
van Marie 2008). The models presented here suggest that the 



most massive stellar winds (M w j n d > 10~ 9 M Q yr -1 ) of young 
T Tauri stars may exist in a similar kind of superposition of 
outflowing and inflowing shells. 

5.3. Analytic Estimate of Wave-Driven Mass Loss 

For the youngest T Tauri models with (seemingly) no steady 
state, it is possible to use an analytic technique to estimate 
how much mass gets accelerated up to the wave-modified crit- 
ical point. As described above, it is not certain whether all of 
this mass can be accelerated to infinity. However, the abil- 
ity to determine a mass flux that applies to the finite region 
between the stellar surface and the critical radius may be suf- 
ficient to predict many observed mass loss diagnostics. Figure 
8a showed that r cr ; t ^> /?* for these young ages, and thus the 
"subcritical volume" is relatively large. 

The modified critical point equation given above for the 
limiting case of dissipationless Alfven waves has been studied 
for several decades (e.g., Jacques 1977; Hartmann & Mac- 
Gregor 1980; DeCampli 1981; Holzer et al. 1983; Wang & 
Sheeley 1991). Equation (|3TT > can be solved for the critical 
point radius by setting the right-hand side to zero, to obtain 



2m 2 ' 

z "crit 



3v 2 



(34) 



where it is assumed that Va (or Ma ^> 1) at the criti- 
cal point, which applies in the youngest T Tauri models. In 
order to compute a value for the critical point radius, how- 
ever, we would have to know the value of v± cr i t , the Alfven 
wave amplitude at the critical point. No straightforward (a pri- 
ori) method was found to predict v^ cr ; t from the other known 
quantities derived in § 3. Thus, we relied on an empirical fit- 
ting relation that was produced from the 6 youngest ZEPHYR 
models (with 0.45 < t < 1.4 Myr), 



V_Lc 



-2.37 



(35) 

18.9 km s" 1 VlOkms-'V 
where v^* is the Alfven wave amplitude at the photosphere. 
The fit is good to 5% accuracy for the 6 youngest numerical 
models, but it is not known whether this level of accuracy 
persists when it is extrapolated to ages younger than ^0.45 
Myr. 

The assumption that the Alfven waves are dissipationless 
allows the conservation of wave action to be used, i.e., 



P v\(u + Va) 2 A 
V a 



constant . 



(36) 



This can be simplified by noting that, for the youngest T Tauri 
models, the photosphere always exhibits u -C Va, and the 
wave-modified critical point always exhibits m > Va. Thus, 
with these assumptions, equation (|36"I > can be applied at these 
two heights to obtain 



/9*Vj_*V A *A» 



PcritVj_ crit M cl - it A crit 



(37) 



VAcrit 

The main unknown is the density p ct i t at the critical point, 
which can be solved in terms of the other known quantities, 

3 / 2 R 2 . /v,.\ 2 

(38) 



Pcrit 
P* 



Vj_crit 



4 7TP*" 2 rit 

and thus the mass flux is determined uniquely at the criti- 
cal point (Holzer et al. 1983). The mass loss rate is M w i n d = 
PcritMcritAcrit, where the critical flux-tube area A C[ j t is normal- 
ized such that as r — > 00, A — > 47rr 2 . 
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Although the above provides a straightforward algorithm 
for estimating M w j„d for a wave-driven wind, there is one other 
unspecified quantity: the isothermal sound speed a at the crit- 
ical point. A completely "cold" wave-driven model would set 
a = 0, but it was found that the extended chromospheres in 
the youngest ZEPHYR models do contribute some gas pres- 
sure to the overall acceleration. The present set of approxi- 
mate models began with an initial guess for the critical point 
temperature (r cr i t w 10 4 K) and then iterated to find a more 
consistent value. The iteration process involved alternating 
between solving for M w i n d at the critical point (as described 
above) and recomputing r cr i t by assuming a balance between 
Alfven-wave turbulent heating and radiative cooling at the 
critical point. The turbulent heating rate Q\ was computed 
as in equation d29l i, and the cooling rate 2iad was assumed 
to remain proportional to its optically thin limit p 2 A{T) (i.e., 
ignoring the ti± correction factor described in § 4.1). The 
temperature at which this balance occurred was found by in- 
verting the tabulated radiative loss function A(T), shown in 
Figure 1 of Cranmer et al. (2007). For each age — between 
13.5 kyr and 0.45 Myr — this process was run for 20 itera- 
tions, but in all cases it converged rapidly within the first 5 
iterations. 

The results of this iterative estimation of M w ind, T CI it, and 
Merit are shown by the dotted curves in Figures 7 and 8. The 
maximum mass loss rate of 2.4 x 10~ 8 M© yr" 1 occurs at an 
age of about 40 kyr. At the youngest ages, the mass loss rate 
declines because the accretion-driven wave power also de- 
clines (see Fig. 5). The most massive wind corresponds to the 
minimum value of the temperature at the critical point, which 
is indistinguishable from the stellar T e fi at that age. At the 
oldest ages for which these solutions were attempted (t > 1 .5 
Myr), the various approximations made above no longer hold 
(e.g., the critical point is no longer super-Alfvenic) and the 
agreement between the analytic estimates and the ZEPHYR 
models breaks down. 

5.4. Varying the Accretion Rate 

As shown in Figure 3 above, the measured values for M acc 
exhibit a wide spread around the mean relation that was used 
to derive the properties of accretion-driven waves. In order 
to explore how variability in the accretion rate affects the 
resulting stellar wind, a set of additional ZEPHYR models 
was constructed with a factor of two decrease and increase 
in M acc in comparison to equation Q. These models were 
constructed mainly for ages around the thermal bifurcation at 
log 1() f 6.25. It was found that the older "hot corona" mod- 
els were insensitive to variations in the accretion rate, even at 
ages where the accretion-driven component to v^* exceeded 
the internal convective component. 

Figure 10 compares the mass loss rate and maximum wind 
temperature for the three sets of ZEPHYR models: the origi- 
nal standard models, ones with half of the accretion rate given 
by equation (01, and ones with double that rate. The specific 
age at which the thermal bifurcation occurs changes by only 
a small amount over the range of modeled accretion rates. 
(Indeed, the standard model and the double-M acc model un- 
dergo thermal bifurcation at nearly the exact same age.) The 
youngest and coolest models show a greater degree of re- 
sponsiveness to the varying accretion rate than do the older 
hot models. The younger models with higher accretion rates 
have larger photospheric MHD wave amplitudes, and thus 
they give rise to larger mass loss rates and cooler extended 
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FIG. 10. — Age-dependent mass loss rates (a) and peak wind tempera- 
tures (b) for ZEPHYR models produced with the standard accretion rate from 
eq. |5J (solid lines), and accretion rates that are double (dashed lines) and half 
(dotted lines) the standard values. Individual ZEPHYR models are denoted 
by symbols. 

chromospheres. 

Figure 1 1 shows how both M acc and M w ; n d vary as a func- 
tion of the photospheric Alfven wave amplitude vj_*. This 
latter quantity is a key intermediary between the mid-latitude 
accretion and the polar mass loss, and thus it is instructive 
to see how both mass flux quantities scale with its evolving 
magnitude. The three sets of ZEPHYR models from Figure 
10 are also shown in Figure 11, and it is noteworthy that there 
is not a simple one-to-one correspondence between M acc and 
Vj_*- The small spread arises because the fundamental stel- 
lar parameters (e.g., and p*) are different at the three ages 
that correspond to the same value of M acc in the three mod- 
els. Thus, it is not surprising that M w j n d is not a "universal" 
function of the wave amplitude either. 

Figure 12 plots the key "wind efficiency" ratio M w ; n d/M aC c 
as a function of age for the three sets of ZEPHYR models 
and the analytic estimates derived in § 5.3. Note that, in con- 
trast to Figure 10, when plotted as a ratio, the models on the 
young/cool side of the thermal bifurcation all seem to collapse 
onto a single curve, while the older/hotter models separate 
(based on differing accretion rates in the denominator). 

For both the youngest ages (f < 0.5 Myr) and for the spe- 
cific case of TW Hya, the model values shown in Figure 12 
seem to agree reasonably well with the observationally de- 
termined ratios that are also shown in Figure 3b. The mod- 
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FIG. 11. — Mass accretion and mass loss rates for the 3 sets of ZEPHYR 
models shown in Figure 10 (with the same line styles) shown versus the pho- 
tospheric Alfven wave amplitude vx*. Mass loss rates for the analytic/cold 
models of § 5.3 are also shown for the standard accretion rate case (dot- 
dashed line). 
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FIG. 12. — Ratio of mass loss rate to mass accretion rate for measured T 
Tauri stars (symbols as in Fig. 3) and for the ZEPHYR and analytic models 
(line styles as in Figs. 10 and 11). 

els having ages between 0.5 and 10 Myr clearly fall well be- 
low the measured mass loss rates. However, even the limited 
agreement with the data is somewhat surprising, since these 
measured values come from the [O I] A6300 forbidden line 
diagnostic that is widely believed to sample the much larger- 
scale disk wind (Hartigan et al. 1995). It is thus possible that 
stellar winds may contribute to observational signatures that 
previously have been assumed to probe only the (disk-related) 
bipolar jets. 

5.5. X-Ray Emission 

Many aspects of the dynamics and energetics of young stars 
and their environments are revealed by high-energy measure- 
ments such as X-ray emission (e.g., Feigelson & Montmerle 
1999). It is thus worthwhile to determine the level of X-ray 
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FIG. 13. — Ratio of X-ray luminosity to total bolometric luminosity for the 
standard run of ZEPHYR models (solid line), and for observations of solar- 
type stars (asterisks), nearby clusters (thin error bars), and the Sun (thick 
error bar); see text for details. 

flux that the modeled polar winds are expected to generate. 
This has been done in an approximate way to produce order- 
of-magnitude estimates, and should be followed up by more 
exact calculations in future work. 

The optically thin radiative loss rate Qrad described in § 4. 1 
was used as a starting point to "count up" the total number of 
photons generated by each radial grid zone in the ZEPHYR 
models. (Finite optical depth effects in Q rdd were ignored here 
because they contribute mainly to temperatures too low to af- 
fect X-ray fluxes.) This rate depends on the plasma tempera- 
ture T, the density p, and the hydrogen ionization fraction x. 
The total radiative loss rate is multiplied by a fraction F that 
gives only those photons that would be observable as X-rays. 
This fraction is estimated for each radial grid zone as 



F = 



Jd\B x (T)S(X) 



(39) 



where B\(T) is the Planck blackbody function and S(X) is 
an X-ray sensitivity function, for which we use that of the 
ROSAT Position Sensitive Proportional Counter (PSPC) in- 
strument, as specified by Judge et al. (2003). The function 
S(X) is nonzero between about 0.1 and 2.4 keV, with a mini- 
mum around 0.3 keV that separates the hard and soft bands. 
The integration in the denominator is taken over all wave- 
lengths. The use of the blackbody function, rather than a true 
optically thin emissivity, was validated by comparison with 
wavelength-limited X-ray radiative loss rates given by Ray- 
mond et al. (1976). Fractions of emissivity (relative to the 
total loss rate) in specific X-ray wavebands were computed 
for T = 0.5, 1, 2, and 5 MK and compared with the plots of 
Raymond et al. (1976). The agreement between the models 
and the published curves was always better than a factor of 
two. 5 

5 Although this is obviously not accurate enough for quantitative compar- 
isons with specific observations, it allows the correct order of magnitude of 
the X-ray emission to be estimated; see Figure 13. Also, the factor-of-two 
validation should be taken in the context of the factor of ~10()0 variation in 
these fractions over the modeled coronal temperatures. 
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Figure 13 shows the simulated ratio of X-ray luminosity Lx 
to the bolometric luminosity L» for the modeled wind regions. 
For each ZEPHYR model, the radiative losses were integrated 
over an assumed spherical volume for the stellar wind (i.e., 
the same assumption used to compute M w ind; see § 5.1) to 
produce Lx- No absorption of X-rays was applied. Figure 13 
also shows a collection of observed X-ray luminosity ratios 
for individual solar-type stars (from Giidel et al. 1998; Garcia - 
Alvarez et al. 2005) and clusters of various ages (Flaccomio 
et al. 2003b; Preibisch & Feigelson 2005; Jeffries et al. 2006). 
For the latter, the error bars indicate the ± la spread about the 
mean values reported in these papers. The range of values for 
the present-day Sun is taken from Judge et al. (2003). 

It is clear from Figure 1 3 that the modeled polar wind re- 
gions do not produce anywhere near enough X-rays to explain 
the observations of young stars. For the present-day Sun, the 
computed value slightly underestimates the lower limit of the 
observed range of X-ray luminosities. This latter prediction 
does make some sense, since the X-ray emission is essen- 
tially computed for the dark "polar coronal holes." On the 
Sun, these never occupy more than about 20% of the sur- 
face. Most notably, though, there is no strong power-law 
decrease in Lx/L* for young ZAMS stars — just as there is 
no power-law decrease in M w j n( j — because there has been no 
attempt to model the rotation-age-activity relationship (Sku- 
manich 1972; Noyes et al. 1984). The ZEPHYR models ne- 
glect the closed-field coronae of these stars (both inside and 
outside the polar-cap regions) that are likely to dominate the 
X-ray emission, like they do for the Sun (e.g., Schwadron et 
al. 2006). 

It is somewhat interesting that the ZEPHYR models un- 
dergo the thermal bifurcation to extended chromospheres (and 
thus disappear from Figure 13 because of the lack of X-ray 
emitting temperatures) at about the same age where there 
seems to exist a slight decline in the X-ray emission of young 
T Tauri stars. Although there is some evidence for such a 
deficit of X-rays in classical T Tauri stars, with respect to the 
older class of weak-lined T Tauri stars (e.g., Flaccomio et al. 
2003a,b; Telleschi et al. 2007), the existence of a distinct peak 
in X-ray emission at intermediate ages has not been proven 
rigorously. If the polar wind regions indeed produce negligi- 
ble X-ray emission compared to the closed-field regions, we 
note that the polar cap area 5 po i grows largest for the youngest 
T Tauri stars (see Fig. 4). Thus, the cooler wind material may 
occupy a significantly larger volume than the hotter (closed- 
loop) coronal plasma at the youngest ages. This may be partly 
responsible for the observed X-ray trends. 

6. DISCUSSION AND CONCLUSIONS 

The primary aim of this paper has been to show how 
accretion-driven waves on the surfaces of T Tauri stars may 
help contribute to the strong rates of atmospheric heating and 
large mass loss rates inferred for these stars. The ZEPHYR 
code, which was originally developed to model the solar 
corona and solar wind (Cranmer et al. 2007), has been ap- 
plied to the T Tauri stellar wind problem. A key aspect of the 
models presented above is that the only true free parameters 
are: (1) the properties of MHD waves injected at the pho- 
tospheric lower boundary, and (2) the background magnetic 
geometry. Everything else (e.g., the radial dependence of the 
rates of chromospheric and coronal heating, the temperature 
structure of the atmosphere, and the wind speeds and mass 
fluxes) emerges naturally from the modeling process. 

For solar-mass T Tauri stars, time-steady stellar winds were 



found to be supportable for all ages older than about 0.45 Myr, 
with accretion rates less than 7 x 10~ 8 M yr" 1 driving mass 
loss rates less than 4 x 10~ 10 Mq yr" 1 . Still younger T Tauri 
stars (i.e., ages between about 13 kyr and 0.45 Myr) may ex- 
hibit time-variable winds with mass loss rates extending up 
several more orders of magnitude to ~ 2 x 10~ 8 M© yr" 1 . 
The transition between time-steady and variable winds occurs 
when the critical point of the flow migrates far enough past 
the Alfven point (at which the wind speed equals the Alfven 
speed) such that the Alfven wave amplitude begins to decline 
rapidly with increasing radius. When this happens, the out- 
ward wave-pressure acceleration is quickly "choked off;" i.e., 
parcels of gas that make it past the critical point cannot be 
accelerated to infinity, and stochastic collisions between up- 
flowing and downflowing parcels must begin to occur. 

The maximum wind efficiency ratio M w i n d/M acc for the T 
Tauri models computed here was approximately 1.4%, com- 
puted for ages of order 0.1 Myr. This is somewhat smaller 
than the values of order 10% required by Matt & Pudritz 
(2005, 2007, 2008) to remove enough angular momentum 
from the young solar system to match present-day conditions. 
It is also well below the observational ratios derived by, e.g., 
Hartigan et al. (1995, 2004), which can reach up to 20% for 
T Tauri stars between 1 and 10 Myr old (see Fig. 12). These 
higher ratios, though, may be the product of both stellar winds 
and disk winds (possibly even dominated by the disk wind 
component). Additionally, it is possible that future observa- 
tional analysis will result in these empirical ratios being re- 
vised upward with more accurate (lower) values of M acc (S. 
Edwards 2008, private communication). 

The accretion-driven solutions for M w ; nc | depend crucially 
on the properties of the waves in the polar regions. It is im- 
portant to note that the calculation of MHD wave properties 
was based on several assumptions that should be examined in 
more detail. The relatively low MHD wave efficiency used 
in equation ( fT7T > is an approximation based on the limiting 
case of waves being far from the impact site. A more realis- 
tic model would have to contain additional information about 
both the nonlinearities of the waves themselves and the verti- 
cal atmospheric structure through which the waves propagate. 
It seems likely that a better treatment of the wave generation 
would lead to larger wave energies at the poles. On the other 
hand, the assumption that the waves do not damp significantly 
between their generation point and the polar wind regions may 
have led to an assumed wave energy that is too high. It is un- 
clear how strong the waves will be in a model that takes ac- 
count of all of the above effects. In any case, the ZEPHYR 
results presented here are the first self-consistent models of 
T Tauri stellar winds that produce wind efficiency ratios that 
even get into the right "ballpark" of the angular momentum 
requirements. 

Additional improvements in the models are needed to make 
further progress. For example, the effects of stellar rota- 
tion should be included, both in the explicit wind dynamics 
(e.g., "magneto-centrifugal driving," as recently applied by 
Holzwarth & Jardine 2007) and in the modified subsurface 
convective activity that is likely to affect photospheric ampli- 
tudes of waves and turbulence. Young and rapidly rotating 
stars are likely to have qualitatively different convective mo- 
tions than are evident in standard (nonrotating, mixing-length) 
models. It is uncertain, though, whether rapid rotation gives 
rise to larger (Kapyla et al. 2007; Brown et al. 2007; Ballot et 
al. 2007) or smaller (Chan 2007) convection eddy velocities 
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at the latitudes of interest for T Tauri stellar winds. In any 
case, rapid rotation can also increase the buoyancy of subsur- 
face magnetic flux elements, leading to a higher rate of flux 
emergence (Holzwarth 2007). Also, the lower gravities of T 
Tauri stars may give rise to a larger fraction of the convective 
velocity reaching the surface as wave energy (e.g., Renzini et 
al. 1977), or the convection may even penetrate directly into 
the photosphere (Montalban et al. 2004). 

Future work must involve not only increased physical re- 
alism for the models, but also quantitative comparisons with 
observations. The methodology outlined in this paper should 
be applied to a set of real stars, rather than to the idealized 
evolutionary sequence of representative parameters. Mea- 
sured stellar masses, radii, and T e ff values, as well as accre- 
tion rates, magnetic field strengths, and hot spot filling factors 
(<5), should be used as constraints on individual calculations 
for the stellar wind properties. It should also be possible to 
use measured three-dimensional magnetic fields (e.g., Donati 
et al. 2007; Jardine et al. 2008) to more accurately map out 
the patterns of accretion stream footpoints, wave fluxes, and 
the flux tubes in which stellar winds are accelerated. 

This work helps to accomplish the larger goal of under- 
standing the physics responsible for low-mass stellar outflows 



and the feedback between accretion disks, winds, and stellar 
magnetic fields. In addition, there are links to more interdisci- 
plinary studies of how stars affect objects in young solar sys- 
tems. For example, the coronal activity and wind of the young 
Sun is likely to have created many of the observed abundance 
and radionuclide patterns in early meteorites (Lee et al. 1998) 
and possibly affected the Earth's atmospheric chemistry to the 
point of influencing the development of life (see, e.g., Ribas 
et al. 2005; Gudel 2007; Cuntz et al. 2008). The identification 
of key physical processes in young stellar winds is important 
not only for understanding stellar and planetary evolution, but 
also for being able to model and predict the present-day im- 
pacts of solar variability and "space weather" (e.g., Feynman 
& Gabriel 2000; Cole 2003). 
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